module PLANAR_PEM
  !
  ! PEM geometry definitions
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
#define green(i) (iand(mana,i)/=0)
  use suty
  use ground_control
  implicit none
  private
  public setenv_ppem
  public printenv_ppem
  public alloc_sysmat
  public set_size_sysmat
  public set_sysmat
  public mSiddon
  public read_list
  public read_hist
  public tros_hist
  public sort_hist
  public load_Cij
  public store_Cij
  public get_minmax_sysmat
  public norm_hist
  public set_hist_mask
  public conv_img
  public median_img
  public resort_sysmat
  public set_idim
  public save_ci
  public load_ci
  public atcUProj
  public flipHistX1,flipHistY1,flipHistB1,flipHistX2

  !
  integer,parameter::VAM_CENTER = 0          ! center to center
  integer,parameter::VAM_EDGE   = 1          ! edge to center
  integer,pointer::voxel_alignment_mode      ! 
  !
  integer,public,parameter::SME_MSIDDON   = 0
  integer,public,parameter::SME_MGAUSS    = 1
  !
  integer,public,parameter::MAT_IMG  = 0
  integer,public,parameter::MAT_PRJ = 1
  !
  integer,parameter,public::OPT_MR_NONE   = 0
  integer,parameter,public::OPT_MR_P1     = 1  
  integer,parameter,public::OPT_MR_P4I    = 4
  integer,parameter,public::OPT_MR_P9I    = 9
  !
  character(MAX_LEN_FILENAME),pointer::fsysmat
  character(MAXPATHLEN),pointer::fhp,fhd,fht
  !
  integer,pointer::opt
  !
  integer,public:: mdxy_max 
  real,parameter::cd=2          ! crystal width in mm
  real,pointer::D               ! distance between crystal surface
  integer,public,pointer::mx    ! number of crystal in x
  integer,public,pointer::my    ! number of crystal in x
  integer,public,pointer::mg    ! number of crystal in a gap
  integer,public,pointer::md    ! number of detector in a head
  integer,public,pointer::mxb   ! number of crystal in x in a detector
  integer,public,pointer::my_true
  integer,public::mxy2
  integer,pointer::mz        ! number of multiray
  integer::ms=1              ! number of sampling points at z
  real,pointer::cdz          ! multiray sampling step in mm
  real,allocatable::cmws(:)
  real,allocatable::cmdx(:),cmdy(:)
  real,allocatable::cwz(:) ! weight in z
  !
  integer,public::nx      ! number of voxels in x = mx*2
  integer,public::ny      ! number of voxels in y = mx*2
  integer,public::ny_true
  integer,public::nz      ! number of voxels in z = arb.
  integer,public::nv      ! nx*ny*nz
  integer,public::nv_true ! nx*ny_true*nz
  integer,public::nh      ! nx/2
  integer,public::nzlo,nzup ! lower upper both positive
  integer jxmin,jxmax,jymin,jymax ! when calculating sysmat
  real,public::dv         ! voxel width
  real,public,pointer::dz ! voxel height in mm
  real,public::jx0,jy0,jz0
  integer,pointer::maxcd
  real max_mdxysqr
  !
  type,public::sysmatel
     integer n         ! number of v
     integer nox,noy
     integer p,q
  end type sysmatel
  !
  integer,parameter,public::SIZE_OF_SYSMATEL=8*5     
  integer,parameter,public::SIZE_OF_SYSMAT_BUF=2*3+4 
  !
  type(sysmatel),allocatable,target,public::sysmat(:,:)
  integer,public::nsysmat=0
  integer,public::nsysmat_allocd=0
  integer,public::bsysmat_allocd=0
  !
  integer mzs
  integer imeth_sysmat
  !
  integer nkern
  real kern_dc
  real,allocatable::kern(:,:,:) ! convolution kernel
  !
  integer,parameter::murphy_factor=3*2
  integer,private::nbuf=0
  !
  real::sysmat_scale=1
  !
contains
  !
  subroutine echo(m,nl)
    character*(*),intent(in)::m
    logical,intent(in),optional::nl
    if(.not.green(OPT_VERBOSE)) return
    if(present(nl)) then
       if(.not.nl) then
          write(*,11) trim(m)
       else
          write(*,12) trim(m)
       end if
    else
       write(*,10) trim(m)
    end if
    call flush(6)
10  format(x,"(pp) ",a)
11  format(x,"(pp) ",a,$)
12  format(a)
  end subroutine echo
  !
  subroutine printenv_ppem
    character*128 str,fig
    if(.not.green(OPT_VERBOSE)) return
    fig=f2a(cd,1)
    str="Detector array: "//trim(i2a(mx))//" x "//trim(i2a(my_true))&
         //" ("//trim(fig)//" x "//trim(fig)//")"
    call echo(str)
    call echo("    Number of blocks in a detector head: "//trim(i2a(md,desc="Not used")))
    call echo("        Number of crystals in a detector head: "//trim(i2a(mxb,desc="Not used")))
    call echo("        Number of crystals in a detector gap: "//trim(i2a(mg,desc="Not used")))
    call echo("Detector separation: "//trim(f2a(D,1)))
    fig=f2a(dv,1)
    str="Voxel: "//trim(i2a(nx))//" x "//trim(i2a(ny_true))//" x "//trim(i2a(nz))&
         //" ("//trim(fig)//" x "//trim(fig)//" x "//trim(f2a(dz,1))//")"
    call echo(str)
    str="FOV: "//trim(f2a(real(nx),1))//" x "//trim(f2a(real(ny_true),1))&
         //" x "//trim(f2a(dz*nz,1))
    call echo(str)
    str="Multiray sampling in z: "//trim(f2a(cdz,1))//" x "//trim(i2a(mz))&
         //" step(s) with "//trim(i2a(ms))//" points."
    call echo(str)
    if(green(OPT_USE_FILTER)) then
       str="Filter decay constant: "//trim(f2a(kern_dc,2))
       call echo(str)
    end if
    str="Voxel alignment:"
    select case(voxel_alignment_mode)
    case(VAM_CENTER)
       str=trim(str)//" Center-to-Center"
    case(VAM_EDGE)
       str=trim(str)//" Edge-to-Center"
    end select
    call echo(str)
    str="Slice range: "//trim(i2a(nzlo))//" -> "//trim(i2a(nzup))
    call echo(str)
  end subroutine printenv_ppem
  !
  subroutine setenv_ppem(pem,w)
    use vuty
    type(ppem_env),intent(in),target::pem
    integer,intent(in)::w
    !
    real,parameter::cmw1p9i=16.0_dp/81.0_dp ! crystal mesh weight
    real,parameter::cmw2p9i=10.0_dp/81.0_dp
    real,parameter::cmw3p9i=25.0_dp/324.0_dp
    !
    real,parameter::cmw1p9b=4.0_dp/9.0_dp
    real,parameter::cmw2p9b=1.0_dp/9.0_dp
    real,parameter::cmw3p9b=1.0_dp/36.0_dp 
    !
    real,parameter::p=sqrt(3.0_dp/5.0_dp)
    real,parameter::q=sqrt(1.0_dp/3.0_dp)
    real,parameter::z=0,u=1
    !
    integer jx,jy,jz
    real rho
    !
    fsysmat=>pem%fsysmat
    fhp=>pem%fhp
    fhd=>pem%fhd
    fht=>pem%fht
    opt=>pem%option
    D=>pem%detector_separation
    maxcd=>pem%maxcd
    if(maxcd>=0) then
       max_mdxysqr=2*maxcd**2!(D*tan(max_incang))**2)/4.d0 !
    else
       max_mdxysqr=-1
    end if
    ! 
    mx=>pem%num_crystal_x
    my=>pem%num_crystal_x ! it's CORRECT
    mdxy_max=min(pem%num_crystal_x,pem%num_crystal_y)
    my_true=>pem%num_crystal_y
    md=>pem%num_detector_in_head
    mg=>pem%num_crystal_in_gap
    mxb=>pem%num_crystal_x_in_detector

    mz=>pem%num_doi_sampling
    if(mz<=0) mz=1
    cdz=>pem%doi_sampling_step
    if(cdz<=0) cdz=1
    dz=>pem%voxel_height
    if(dz<=0) call echo_red("*** Invalid voxel height: "//trim(f2a(dz,-1)))
    !
    ! w must be 2 or 1.
    voxel_alignment_mode=>pem%vam
    select case(voxel_alignment_mode)
    case(VAM_CENTER)
       nx=mx*w+1
       ny=my*w+1
       ny_true=my_true*w+1
       jxmin=1
       jymin=1
       jxmax=nx
       jymax=ny
    case(VAM_EDGE)
       nx=mx*w     
       ny=my*w     
       ny_true=my_true*w
       jxmin=1
       jymin=1
       jxmax=nx
       jymax=ny
    case default
       call echo_red("*** Invalid voxel alignment")
    end select
    !
    nh=max(nx,ny)
    select case(w)
    case(2)
       dv=cd/2.0_dp
       if(iand(nh,1)==1) then
          nh=nh/2+1
       else
          nh=nh/2
       end if
    case(1)
       dv=cd
    case default
       call echo_red("*** Voxel width must be 1 or 1/2 of crystal.")
    end select
    nz=int(D/dz-1)   ! <<<<<<<<<<<< for cubic FOV to rotate
    if(iand(nz,1)==0) nz=nz-1
    nv=nx*ny*nz
    nv_true=nx*ny_true*nz
    mxy2=(mx*my)**2
    !
    if(pem%lo==-huge(0)) then
       nzlo=1
    else
       nzlo=max(0,pem%lo+nz/2)
    end if
    if(pem%up==huge(0)) then
       nzup=nz
    else
       nzup=min(nz,pem%up+nz/2)
    end if
    !
    ! x = (jx - jx0) * dv
    jx0=real(nx+1)/2.0_dp
    jy0=real(ny+1)/2.0_dp
    jz0=real(nz+1)/2.0_dp
    !
    call init_cmbuf
    select case(pem%multiray_option)
    case(OPT_MR_P9I)
       ms=9
       call alloc_cmbuf
       cmws=[cmw1p9i,cmw2p9i,cmw2p9i,cmw2p9i,cmw2p9i,cmw3p9i,&
            cmw3p9i,cmw3p9i,cmw3p9i]
       cmdx=[z, p, z, -p,  z, p, -p, -p,  p]
       cmdy=[z, z, p,  z, -p, p,  p, -p, -p]
    case(OPT_MR_P4I)
       ms=4
       call alloc_cmbuf
       cmws=1.0_dp/4.0_dp
       cmdx=[z, q, -q, -q,  q]
       cmdy=[z, q,  q, -q, -q]
    case(OPT_MR_NONE)
    case(OPT_MR_P1)
    case default
       call echo_red("*** Invalid multiray_option.")
    end select
    cmws=cmws/mz
    !
    mzs=mz*ms
    !
    if(allocated(cwz)) deallocate(cwz)
    allocate(cwz(nz))
    do jz=1,nz
       cwz(jz)=get_cwz((real(jz)-jz0)*dz)
    end do
    !
    select case(pem%sysmat_evaluator)
    case(SME_MSIDDON,SME_MGAUSS)
       imeth_sysmat=pem%sysmat_evaluator
    case default
       call echo_red("*** Invalid sysmat_evaluator.")
    end select
    !
    if(allocated(kern)) deallocate(kern)
    nkern=0
    if(pem%fa>0) then
       nkern=1
       kern_dc=4*log(2.0_dp)/pem%fa**2
       allocate(kern(-nkern:nkern,-nkern:nkern,-nkern:nkern))
       do jz=-nkern,nkern
          do jy=-nkern,nkern
             do jx=-nkern,nkern
                rho=sqrt(jx**2.0_dp+jy**2.0_dp+(jz*dz)**2.0_dp)
                kern(jx,jy,jz)=exp(-kern_dc*rho)
             end do
          end do
       end do
    end if
    !
    call initenvv(mx,nh)
    !
  contains
    !
    real function get_cwz(z)
      real,intent(in)::z      
      ! z in mm
      get_cwz=1
!      get_cwz=-0.03*(z/10.d0)**2-0.03*abs(z/10.d0)+1
    end function get_cwz
    !
    subroutine init_cmbuf
      ms=1
      call alloc_cmbuf
      cmws=1
      cmdx=0
      cmdy=0
    end subroutine init_cmbuf
    !
    subroutine dealloc_cmbuf
      if(allocated(cmws)) deallocate(cmws)
      if(allocated(cmdx)) deallocate(cmdx)
      if(allocated(cmdy)) deallocate(cmdy)
    end subroutine dealloc_cmbuf
    !
    subroutine alloc_cmbuf
      call dealloc_cmbuf
      allocate(cmws(ms))
      allocate(cmdx(ms),cmdy(ms))
    end subroutine alloc_cmbuf
  end subroutine setenv_ppem
  !
  subroutine set_LOR_ends(ix,iy,X,w)
    integer,intent(in)::ix,iy
    real,intent(out)::X(3,ms*mz)
    real,intent(out)::w(ms*mz)
    integer iz,k,j
    real X0(3)
    X0(1)=(cd*ix-mx)-1
    X0(2)=(cd*iy-my)-1
    do iz=1,mz
       X0(3)=D/2.0_dp+(iz-1)*cdz
       do k=1,ms
          j=(iz-1)*ms+k
          X(1,j)=X0(1)+cmdx(k)
          X(2,j)=X0(2)+cmdy(k)
          X(3,j)=X0(3)
          w(j)=cmws(k)
       end do
    end do
  end subroutine set_LOR_ends
  !
  subroutine alloc_sysmat_buf(n,sm)
    integer,intent(in)::n
    type(sysmatel),intent(out)::sm
    sm%n=n
    sm%p=malloc(n*SIZE_OF_SYSMATEL)
  end subroutine alloc_sysmat_buf
  !
  subroutine echo_red(mess,code)
    character*(*),intent(in)::mess
    integer,intent(in),optional::code
    integer ec
    character*(*),parameter::swan_song="*** emergency escape from module `planar_pem'."
    write(0,*) trim(adjustl(mess))
    call set_rstat(RSTAT_ERROR)
    ec=-1
    if(present(code)) ec=code
    call exit(ec)
  end subroutine echo_red
  !
  subroutine save_hist(f,h)
    character*(*),intent(inout)::f
    integer*2,intent(in),target::h(mx,my,mx,my)
    integer*2,pointer::hout(:,:,:,:)
    integer iu,istat
    integer ix1,iy1,ix2,iy2
    iu=proc_output(f,20,istat,".dat")
    if(iu>0) then
       if(my/=my_true) then
          allocate(hout(mx,my_true,mx,my_true))
          hout=0_2
          forall(ix1=1:mx,iy1=1:my_true,ix2=1:mx,iy2=1:my_true)
             hout(ix1,iy1,ix2,iy2)=h(ix1,iy1,ix2,iy2)
          end forall
       else
          hout=>h
       end if
       call echo("Writing histogram: "//trim(f),.false.)
       write(iu,iostat=istat) hout
       close(iu)
       if(my/=my_true) deallocate(hout)
       call echo(" ["//trim(iostatus(istat))//"]",.true.)
    else
       call echo_red(mess_ioerr(EID_OPEN,f,istat))
    end if
  end subroutine save_hist    
  !
  subroutine read_list(n,hist,indx,file)
    use buty
    integer,intent(in)::n
    real*4,intent(out)::hist(mx,my,mx,my)
    integer,intent(in),optional::indx(n)
    integer*2,target::whistpd(mx,my,mx,my,2)
    integer*2,pointer::whp(:,:,:,:),whd(:,:,:,:)
    character*(*),intent(in),optional::file    
    integer plist,szlist
    integer*8 qw
    character*128 mess
    pointer(p,qw)
    call mkhis_init(my,mxb,mg,md,n,FT_LIST,green(OPT_VERBOSE),indx,file)
    if(shmstats(SHMID_IN)%p>0) then
       p=shmstats(SHMID_IN)%p
       szlist=qw ! file size
       p=p+8     ! data
       plist=p
    else
       plist=0
       szlist=0
    end if
    whistpd=0_2
    whp=>whistpd(:,:,:,:,2)
    whd=>whistpd(:,:,:,:,1)
    call mkhisqh(plist,szlist,loc(whistpd))
    if(associated(fhp)) call save_hist(fhp,whp)
    if(associated(fhd)) call save_hist(fhd,whd)
    if(.not.green(OPT_USE_TRUE)) then
       hist=whp
       mess="prompt coin."
    end if
    if(green(OPT_USE_TRUE).or.associated(fht)) then
       whp=whp-whd
       where(whp<0) whp=0_2
       if(associated(fht)) call save_hist(fht,whp)
       if(green(OPT_USE_TRUE)) then
          hist=whp       
          mess="delayed subtraced."
       end if
    end if
    mess="Histogram loaded with "//trim(mess)
    call echo(mess)
  end subroutine read_list
  !
  subroutine load_hist(hist)
    real*4,intent(out)::hist(mx,my,mx,my)
    integer*2,allocatable::buf2(:,:,:,:)
    integer*4,allocatable::buf4(:,:,:,:)
    if(green(OPT_READ_32)) then
       allocate(buf4(mx,my_true,mx,my_true))
       call mcps(loc(buf4),shmstats(SHMID_IN)%p,mx*my_true*mx*my_true*4)
       call set_hist_32(buf4,hist)
    else
       allocate(buf2(mx,my_true,mx,my_true))
       call mcps(loc(buf2),shmstats(SHMID_IN)%p,mx*my_true*mx*my_true*2)
       call set_hist_16(buf2,hist)
    end if
    where(hist<0) hist=0_4
  end subroutine load_hist
  !
  subroutine read_hist(n,hist,indx,file)
    integer,intent(in)::n
    real*4,intent(out)::hist(mx,my,mx,my)
    integer,intent(in),optional::indx(n)
    character*(*),intent(in),optional::file
    character(MAXPATHLEN) fhist
    integer*2,allocatable::b2(:,:,:,:)
    integer*4,allocatable::b4(:,:,:,:)
    integer*2,allocatable::buf2(:,:,:,:)
    integer*4,allocatable::buf4(:,:,:,:)
    integer istat,ifile
    integer nf
    if(green(OPT_READ_32)) then
       allocate(buf4(mx,my_true,mx,my_true))       
       allocate(b4(mx,my_true,mx,my_true))
       allocate(b2(0,0,0,0))
       b4=0_4
    else
       allocate(buf2(mx,my_true,mx,my_true)) 
       allocate(b2(mx,my_true,mx,my_true))
       allocate(b4(0,0,0,0))
       b2=0_2
    end if
    nf=n
    if(present(file)) nf=1
    do ifile=1,nf
       if(present(indx)) then
          if(indx(ifile)/=FT_HIST) cycle
          call get_command_argument(ifile,fhist)
       else if(present(file)) then
          fhist=trim(file)
       end if
       open(unit=1,file=fhist,access="stream",iostat=istat,status="old")
       if(istat/=0) call echo_red(mess_ioerr(EID_OPEN,fhist,istat))
       if(green(OPT_READ_32)) then
          read(1,iostat=istat) buf4
          b4=b4+buf4
       else
          read(1,iostat=istat) buf2
          b2=b2+buf2
       end if
       if(istat/=0) call echo_red(mess_ioerr(EID_READ,fhist,istat))         
       call echo("Histogram loaded: "//trim(fhist))
       close(1)
    end do
    if(green(OPT_READ_32)) then
       call set_hist_32(b4,hist)
       deallocate(b4,buf4,b2)
    else
       call set_hist_16(b2,hist)
       deallocate(b2,buf2,b4)
    end if
    where(hist<0.0_4) hist=0.0_4
  end subroutine read_hist
  !
  subroutine set_hist_32(b,h)
    integer*4,intent(in)::b(mx,my_true,mx,my_true)
    real*4,intent(out)::h(mx,my,mx,my)
    integer ix1,iy1,ix2,iy2
    if(my==my_true) then
       h=b
    else
       h=0.0_4
       forall(ix1=1:mx,iy1=1:my_true,ix2=1:mx,iy2=1:my_true)
          h(ix1,iy1,ix2,iy2)=b(ix1,iy1,ix2,iy2)
       end forall
    end if    
  end subroutine set_hist_32
  !
  subroutine set_hist_16(b,h)
    integer*2,intent(in)::b(mx,my_true,mx,my_true)
    real*4,intent(out)::h(mx,my,mx,my)    
    integer ix1,iy1,ix2,iy2
    if(my==my_true) then
       h=b
    else
       h=0.0_4
       forall(ix1=1:mx,iy1=1:my_true,ix2=1:mx,iy2=1:my_true)
          h(ix1,iy1,ix2,iy2)=b(ix1,iy1,ix2,iy2)
       end forall       
    end if
  end subroutine set_hist_16
  !
  subroutine tros_hist(h)
    real*4,intent(inout)::h(:,:,:,:)
    real*4,allocatable::buf(:,:,:,:)
    integer ix1,iy1,ix2,iy2
    integer dx,dy
    allocate(buf(mx,my,mx,my))
    do iy2=1,my
       do ix2=1,mx
          do iy1=1,my
             dy=iy2-iy1
             if(dy>=0) then
                dy=dy+1
             else
                dy=my+dy+1
             end if
             do ix1=1,mx
                dx=ix2-ix1
                if(dx>=0) then
                   dx=dx+1
                else
                   dx=mx+dx+1
                end if
                buf(ix1,iy1,dx,dy)=h(ix1,iy1,ix2,iy2)
             end do
          end do
       end do
    end do
    h=buf
    deallocate(buf)
  end subroutine tros_hist
  !
  subroutine sort_hist(h)
    real*4,intent(inout)::h(:,:,:,:)
    real*4,allocatable::buf(:,:,:,:)
    integer ix1,iy1,ix2,iy2
    integer dx,dy
    allocate(buf(mx,my,mx,my))
    do iy2=1,my
       do ix2=1,mx
          do iy1=1,my
             dy=iy2-iy1
             if(dy>=0) then
                dy=dy+1
             else
                dy=my+dy+1
             end if
             do ix1=1,mx
                dx=ix2-ix1
                if(dx>=0) then
                   dx=dx+1
                else
                   dx=mx+dx+1
                end if
                buf(ix1,iy1,ix2,iy2)=h(ix1,iy1,dx,dy)
             end do
          end do
       end do
    end do
    h=buf
    deallocate(buf)
  end subroutine sort_hist
  !
  subroutine norm_hist(ih,oh,nh,nf)
    real*4,intent(in)::ih(:,:,:,:)
    real*4,intent(out)::oh(:,:,:,:)
    real*4,intent(in),optional,target::nh(:,:,:,:)
    character*(*),intent(in),optional::nf
    integer*2,allocatable::buf(:,:,:,:)
    integer*4,allocatable::buf4(:,:,:,:)
    real*4,pointer::norm(:,:,:,:)
    integer ix1,iy1,ix2,iy2,n,istat
    real s
    if(present(nf)) then
       open(unit=1,file=nf,access="stream",iostat=istat,status="old")
       if(istat/=0) call echo_red(mess_ioerr(EID_OPEN,nf,istat))
       if(green(OPT_READ_32)) then
          allocate(buf4(mx,my,mx,my))
          read(1,iostat=istat) buf4
       else
          allocate(buf(mx,my,mx,my))
          read(1,iostat=istat) buf
       end if
       close(1)
       if(istat/=0) call echo_red(mess_ioerr(EID_READ,nf,istat))
       call echo("Normalize file loaded: "//trim(nf))
       allocate(norm(mx,my,mx,my))
       if(green(OPT_READ_32)) then
          norm=buf4
          deallocate(buf4)
       else
          norm=buf
          deallocate(buf)
       end if
       call tros_hist(norm)
    else if(present(nh)) then
       norm=>nh
    else
       return
    end if
    s=0
    n=0
    do iy2=1,my
       do ix2=1,mx
          do iy1=1,my
             do ix1=1,mx
                if(norm(ix1,iy1,ix2,iy2)==0) cycle
                s=s+norm(ix1,iy1,ix2,iy2)
                n=n+1
             end do
          end do
       end do
    end do
    s=s/n
    norm=real(norm/s,4)
    where(norm==0)
       oh=0.0_4
    elsewhere
       oh=ih/norm
    end where
    if(present(nf)) deallocate(norm)
  end subroutine norm_hist
  !
  subroutine set_hist_mask(oh)
    real*4,intent(inout)::oh(:,:,:,:)
    integer id,ig,ix
    do id=1,md-1
       do ig=1,mg
          ix=id*mxb+mg*(id-1)+ig
          oh(ix,1:my,1:mx,1:my)=0.0_4
          oh(1:mx,1:my,ix,1:my)=0.0_4
       end do
    end do

  end subroutine set_hist_mask
  !
  integer function load_Cij(fin)
    character*(*),intent(in)::fin
    integer ix,iy,n,istat
    real D_,dz_
    integer mx_,my_,nx_,ny_,nz_
    integer*2 jxs(*),jys(*),jzs(*)
    real*4 ws(*)
    pointer(pjxs,jxs)
    pointer(pjys,jys)
    pointer(pjzs,jzs)
    pointer(pws,ws)
    !
    call alloc_sysmat
    call echo( "Opening "//trim(fin)//"...")
    open(unit=1,file=fin,access="stream",iostat=istat,status="old")
    if(istat/=0) call echo_red(mess_ioerr(EID_OPEN,fin,istat))
    read(1,iostat=istat) D_,mx_,my_,nx_,ny_,nz_,dz_
    if(istat/=0) call echo_red(mess_ioerr(EID_READ,fin,istat))
    if(D_==D.and.mx_==mx.and.my_==my.and.nx_==nx.and.ny_==ny.and.nz_==nz.and.dz_==dz) then
       do iy=1,my
          do ix=iy,mx
             read(1,iostat=istat) sysmat(ix,iy)%n
          end do
       end do
       if(istat/=0) call echo_red(mess_ioerr(EID_READ,fin,istat))
       do iy=1,my
          do ix=iy,mx
             n=sysmat(ix,iy)%n
             if(n==0) cycle
             call alloc_sysmat_buf(n,sysmat(ix,iy))
             pjxs=sysmat(ix,iy)%p
             pjys=pjxs+n*2
             pjzs=pjys+n*2
             pws =pjzs+n*2
             read(1,iostat=istat) jxs(1:n),jys(1:n),jzs(1:n),ws(1:n)
             if(istat/=0) exit
          end do
       end do
       if(istat/=0) call echo_red(mess_ioerr(EID_READ,fin,istat))
       load_Cij=0
    else
       write(0,*) "*** Error geometry mismatch detected:"
       write(0,100) D_,D,mx_,mx,my_,my,nx_,nx,ny_,ny,nz_,nz,dz_,dz
       100 format(5x,"  =source=   =target=",/,&
                x,"D ",4x,f8.4,3x,f8.4,/,&
                x,"mx",9x,i3,8x,i3,/,&
                x,"my",9x,i3,8x,i3,/,&
                x,"nx",9x,i3,8x,i3,/,&
                x,"ny",9x,i3,8x,i3,/,&
                x,"nz",9x,i3,8x,i3,/,&
                x,"dz",4x,f8.4,3x,f8.4)
       load_Cij=1
    end if
    close(1)
    call echo("read_Cij returns with code = "//trim(i2a(load_Cij)))
  end function load_Cij
  !
  subroutine store_Cij(fout)
    character*(*),intent(in)::fout
    character*256 f
    integer ix,iy,n
    integer iu,istat
    integer*2 jxs(*),jys(*),jzs(*)
    real*4 ws(*)
    pointer(pjxs,jxs)
    pointer(pjys,jys)
    pointer(pjzs,jzs)
    pointer(pws,ws)
    !
    f=trim(fout)//".dat"
    call echo("Writing "//trim(adjustl(f)))
    iu=proc_output(f,20,istat)
    if(iu==-1) call echo_red(mess_ioerr(EID_OPEN,f,istat))
    write(iu,iostat=istat) D,mx,my,nx,ny,nz,dz
    if(istat/=0) call echo_red(mess_ioerr(EID_WRITE,f,istat))
    do iy=1,my
       do ix=iy,mx
          write(iu,iostat=istat) sysmat(ix,iy)%n
          if(istat/=0) exit
       end do
    end do
    if(istat/=0) call echo_red(mess_ioerr(EID_WRITE,f,istat))
    do iy=1,my
       do ix=iy,mx
          n=sysmat(ix,iy)%n
          if(n==0) cycle
          pjxs=sysmat(ix,iy)%p
          pjys=pjxs+n*2
          pjzs=pjys+n*2
          pws =pjzs+n*2
          write(iu,iostat=istat) jxs(1:n),jys(1:n),jzs(1:n),ws(1:n)
          if(istat/=0) exit
       end do
    end do
    close(iu)
    if(istat/=0) call echo_red(mess_ioerr(EID_WRITE,f,istat))
  end subroutine store_Cij
  !
  subroutine dealloc_sysmat
    integer ix,iy
    do iy=1,my
       do ix=iy,mx
          if(sysmat(ix,iy)%p/=0) then
             call free(sysmat(ix,iy)%p)
             sysmat(ix,iy)%p=0
          end if
          if(sysmat(ix,iy)%q/=0) then
             call free(sysmat(ix,iy)%q)
             sysmat(ix,iy)%q=0
          end if
       end do
    end do
    deallocate(sysmat)
  end subroutine dealloc_sysmat
  !
  subroutine init_sysmat
    sysmat%n=0
    sysmat%nox=0
    sysmat%noy=0
    sysmat%p=0
    sysmat%q=0
  end subroutine init_sysmat
  !
  subroutine alloc_sysmat
    if(allocated(sysmat)) call dealloc_sysmat
    allocate(sysmat(mx,my))
    call init_sysmat
    nsysmat_allocd=0
    call echo("Size of raw sysmat: "//trim(i2a(sizeof(sysmat))))
  end subroutine alloc_sysmat
  !
  ! Fails with bus error when optimized
!!$  subroutine scale_sysmat(norm)
!!$    real,intent(in)::norm
!!$    integer ix,iy,k
!!$    integer n
!!$    real*4 w(*)
!!$    real*4 f
!!$    pointer(pw,w)
!!$    f=norm
!!$    do iy=1,my
!!$       do ix=iy,mx
!!$          n=sysmat(ix,iy)%n
!!$          if(n==0) cycle
!!$          pw=sysmat(ix,iy)%p+n*2*3
!!$          w(1:n)=w(1:n)*f
!!$       end do
!!$    end do
!!$    call echo("sysmat scaled by: "//trim(f2a(norm,-1)))
!!$  end subroutine scale_sysmat
  !
  subroutine get_minmax_sysmat(wmin,wmax)
    real,intent(out)::wmin,wmax
    integer ix,iy
    integer k,n
    real*4 w(*)
    pointer(pw,w)
    wmin=huge(1.0_dp)
    wmax=-huge(1.0_dp)
    do iy=1,my
       do ix=iy,mx
          n=sysmat(ix,iy)%n
          pw=sysmat(ix,iy)%p+n*2*3
          do k=1,n
             if(w(k)>wmax)then
                wmax=w(k)
             else if(w(k)<wmin) then
                wmin=w(k)
             end if
          end do
       end do
    end do
  end subroutine get_minmax_sysmat
  !
  subroutine set_size_sysmat
    integer ix,iy
    do iy=1,my
       do ix=iy,mx
          nsysmat_allocd=nsysmat_allocd+sysmat(ix,iy)%n
       end do
    end do
    bsysmat_allocd=nsysmat_allocd*SIZE_OF_SYSMAT_BUF
    call echo("Allocated number of sysmat element: "//trim(i2a(nsysmat_allocd)))
    call echo("Size of allocated matrix (Byte): "//trim(i2a(bsysmat_allocd)))
    call echo("Maximum number of sub-element: "//trim(i2a(maxval(sysmat%n))))
    call echo("Scaling factor of sysmat: "//trim(f2a(sysmat_scale,-1)))
  end subroutine set_size_sysmat
  !
  integer function mSiddon(dy2)
    integer,intent(in)::dy2
    integer dx2
    real,pointer::X1(:)   ! start of LOR
    real,pointer::X2(:)   ! end of LOR
    real d12     ! length of LOR
    real dX12(3) ! X2-X1
    real Xp1(3)  ! coordinate of the first plane 
    real XpN(3)  ! coordinate of the last plane 
    real dX1p1(3)
    real rXp1_1(3)
    ! a must be larger than 0 and less than 1.
    real::a1(3)=0   ! first intersection of LOR in x and y
    real:: aN(3)=0  ! last intersection of LOR in x and y 
    real amin       ! first intersection of LOR
    real amax       ! last intersection of LOR
    integer imin(3) ! first index of plane
    integer imax(3) ! last index of plane
    real,allocatable::a(:)!(0:nabuf)    ! intersection
    real wbuf(nbuf),wwbuf(nbuf)
    integer jbuf(nbuf),per(nbuf)
    integer Ns(3)
    ! etc
    integer i,k,m,n,ic
    integer,target::js(3)
    integer,pointer::jx,jy,jz
    real l,am
    logical skip(3)
    integer ix2,iy2
    real::aa1(3)=0,aa2(3)=0
    real mins(3),maxs(3)
    integer kbuf
    real norm
    real,target::X1s(3,mzs),X2s(3,mzs)
    real,target::w1s(mzs),w2s(mzs)
    integer izs1,izs2
    integer*2 ax
    real ds(3)
    pointer(pax,ax)
    real,parameter::wmin=1.e-9_dp
    !--------------------------------!
    real*4 ws(*)
    real w
    integer*2 jxs(*),jys(*),jzs(*)
    pointer(pjx,jxs)
    pointer(pjy,jys)
    pointer(pjz,jzs)
    pointer(pw,ws)
    !--------------------------------!
    logical fake
    integer nabuf
    !
    nabuf=max(nx+2,ny+2,nz+2)*3
    allocate(a(0:nabuf))
    fake=(nbuf==0)
    mSiddon=0
    call init
    iy2=1+dy2
    main: do dx2=my_true-1,dy2,-1
       ix2=1+dx2
       call set_LOR_ends(ix2,iy2,X2s,w2s)
       kbuf=0
       do izs2=1,mzs
          X2=>X2s(:,izs2)
          do izs1=1,mzs
             X1=>X1s(:,izs1)
             d12=sqrt(sum((X2-X1)**2))
             dX12=X2-X1
             dX1p1=X1-Xp1
             rXp1_1=(Xp1-X1)/dX12
             norm=d12
             where(dX12/=0)
                a1=(Xp1-X1)/dX12
                aN=(XpN-X1)/dX12
                skip=.false.
             elsewhere
                a1=-1
                aN=-1
                skip=.true.
             end where
             maxs=max(a1,aN)
             mins=min(a1,aN)
             where(maxs>1) maxs=0
             where(mins<0) mins=1
             amin=minval(mins)
             amax=maxval(maxs)
             where(dX12>0)
                aa1=amin
                aa2=amax
             elsewhere
                aa1=amax
                aa2=amin
             end where
             imin=int(Ns-(XpN-aa1*dX12-X1)/ds)
             imax=int(1+(X1+aa2*dX12-Xp1)/ds) 
             where(skip)
                imin=0
                imax=0
             end where
             imin=max(1,imin)
             imax=min(Ns,imax)
             n=sum(imax-imin)+3
             k=0
             a(0)=amin
             do ic=1,3
                if(skip(ic)) cycle
                do i=imin(ic),imax(ic)
                   k=k+1
                   if(k>nabuf) then
                      mSiddon=EID_BUFOVER_1
                      exit main
                   end if
                   a(k)=rXp1_1(ic)+(i-1)/dX12(ic)*ds(ic)
                   if(a(k)>=amax) then
                      k=k-1
                      exit
                   elseif(a(k)<=amin) then
                      k=k-1
                   end if
                end do
             end do
             k=k+1
             a(k)=amax
             n=k
             call merge_a
             if(n<1) exit
             !
             k=kbuf
             do m=1,n
                l=a(m)-a(m-1)
                am=(a(m)+a(m-1))/2.0_dp
                js=int(1.0_dp+(dX1p1+am*dX12)/ds )
                w=l*norm*cwz(jz)
                IF(STORE_JSW()==-1) EXIT MAIN
             end do
             kbuf=k
          end do
       end do
       n=kbuf
       if(fake) exit main
       if(n>0) then
          call add_a
          call alloc_sysmat_buf(n,sysmat(ix2,iy2))
          call unpack2(n,loc(jbuf),sysmat(ix2,iy2)%p,loc(wbuf))
       end if
    end do main
    !
    deallocate(a)
    if(fake) mSiddon=n*murphy_factor
    !
  contains 

    integer function store_jsw()
      store_jsw=0
      if(jz<nzlo.or.jz>nzup &
           .or.jx<jxmin.or.jx>jxmax&
           .or.jy<jymin.or.jy>jymax) return
      if(w<=wmin) return
      k=k+1
      if(fake) return
      if(k>nbuf) then
         mSiddon=EID_BUFOVER_2
         STORE_JSW=-1
         return
      end if
      IF(K>NBUF) call echo_red("*** NBUF OVERFLOW!!!")
      pax=loc(jbuf(k))
      ax=transfer(jx,ax)
      pax=pax+2
      ax=transfer(jy,ax)
      pax=pax+2
      ax=transfer(jz,ax)
      wbuf(k)=w*sysmat_scale
    end function store_jsw
    !
    subroutine init
      Xp1(1)=-real(nx)*dv/2.0_dp
      Xp1(2)=-real(ny)*dv/2.0_dp
      Xp1(3)=-real(nz)*dz/2.0_dp
      XpN=-Xp1
      !
      jx=>js(1)
      jy=>js(2)
      jz=>js(3)
      Ns(1)=nx+1
      Ns(2)=ny+1
      Ns(3)=nz+1
      !
      call set_LOR_ends(1,1,X1s,w1s)
      X1s(3,:)=-X1s(3,:)
      !
      ds(1:2)=dv
      ds(3)  =dz
    end subroutine init
    !
    subroutine merge_a
      integer src,dst
      call quickSort_r8(a(1:n),1,n)
      src=0
      dst=0
      do while(src<n)
         src=src+1
         if(a(src)/=a(dst)) then
            dst=dst+1
            if(dst/=src) a(dst)=a(src)
         end if
      end do
      n=dst
    end subroutine merge_a
    !
    subroutine add_a
      integer src,trg,dst
      call persort_i(n,jbuf(1:n),per(1:n))
      src=1
      trg=1
      dst=1
      wwbuf(1:n)=wbuf(1:n)
      do while(src<n)
         src=src+1
         if(jbuf(src)/=jbuf(trg)) then
            dst=dst+1
            trg=src
            wwbuf(dst)=wbuf(src)
            jbuf(dst)=jbuf(src)
         else
            wwbuf(dst)=wwbuf(dst)+wbuf(per(src))
         end if
      end do
      n=dst
      wbuf(1:n)=wwbuf(1:n)
    end subroutine add_a
    !
  end function mSiddon
  !
  subroutine set_sysmat(scale)
    use omp_lib
    interface 
       integer function eval_sysmat(a)
         integer,intent(in)::a
       end function eval_sysmat
    end interface    
    real,intent(in),optional::scale
    integer dy2
    real t0,t1
    integer rts(8,2)
    integer rc
    real wmin,wmax
    character*32 str
    integer ims
    pointer(ptr_eval_sysmat,eval_sysmat)
    !
    if(present(scale)) then
       sysmat_scale=scale
       ims=imeth_sysmat
    else
       sysmat_scale=1
       ims=SME_MSIDDON
    end if
    !
    if(green(OPT_LOAD_SYSMAT)) then
       rc=load_Cij(fsysmat)
       if(rc==0.and.green(OPT_SAVE_SYSMAT)) &
            call gc_veto_C(OPT_SAVE_SYSMAT)
    else
       rc=-1
    end if
    !
    if(rc/=0) then
       str="MR Siddon"
       ptr_eval_sysmat=loc(mSiddon)
       select case(ims)
       case(SME_MGAUSS)
          str="MR Siddon G"
          ptr_eval_sysmat=loc(mSiddong)
       end select
       call echo("Sysmat evaluator: "//trim(str))
       nbuf=0
       nbuf=eval_sysmat(my_true-1)
       t0=time8()
       call date_and_time(values=rts(:,1))
       call alloc_sysmat
       !$omp parallel do schedule(dynamic) private(dy2)
       do dy2=0,my_true-1
          select case(eval_sysmat(dy2))
          case(EID_BUFOVER_1)
             call echo_red("*** INTERNAL ERROR bufferoverflow a: "&
                  //trim(i2a(nbuf)))
          case(EID_BUFOVER_2)
             call echo_red( "*** INTERNAL ERROR bufferoverflow b: "&
                  //trim(i2a(nbuf)))
          end select
       end do
       !$omp end parallel do
       call date_and_time(values=rts(:,2))
       t1=time8()
       call echo("set_sysmat elapsed time: "//trim(f2a(t1-t0,3))//"u "&
            //trim(f2a(etime(rts(:,1),rts(:,2)),3))//"s")
    end if
    call set_size_sysmat
    call get_minmax_sysmat(wmin,wmax)
    call echo("sysmat min/max: "//trim(f2a(wmin,-1))//" / "//trim(f2a(wmax,-1)))
    if(green(OPT_SAVE_SYSMAT)) call store_Cij(fsysmat)
  end subroutine set_sysmat
  !
  subroutine median_img(img,m)
    real*4,intent(inout)::img(:,:,:)
    integer,intent(in)::m
    integer n    
    real*4 buf(nx,ny,nz)
    real*4 v(-m:m,-m:m)
    real*4 vv(*)
    pointer(p,vv)
    integer jx,jy,jz
    integer kx,ky
    buf=0.0_sp
    p=loc(v)
    n=(2*m+1)**2
    do jz=1,nz
       do jy=1+m,ny-m
          do jx=1+m,nx-m             
             do ky=-m,m
                do kx=-m,m
                   v(kx,ky)=img(jx+kx,jy+ky,jz)
                end do
             end do
             buf(jx,jy,jz)=median(n,vv)
          end do
       end do
    end do
    img=buf
  end subroutine median_img
  !
  subroutine conv_img(img,is3d)
    real*4,intent(inout)::img(:,:,:)
    logical,intent(in)::is3d
    real*4 tmpimg(nx,ny,nz)
    real sin,sout
    integer kz1,kz2
    integer jz_
    if(nkern==0) return
    if(is3d) then
       kz1=-nkern
       kz2=nkern
    else
       kz1=0
       kz2=0
    end if
    sin=sum(img)
    if(sin==0) return
    tmpimg=0.0_4
    !$omp parallel do private(jz_) schedule(dynamic)
    do jz_=1,nz
       call do_conv(jz_)
    end do
    !$omp end parallel do
    sout=sum(tmpimg)
    if(sout==0) return
    img=tmpimg*real(sin/sout,4)
  contains
    subroutine do_conv(jz)
      integer,intent(in)::jz
      integer jx,jy
      integer kx,ky,kz
      integer jjx,jjy,jjz
      real s
      do jy=1,ny
         do jx=1,nx
            s=0
            do kz=kz1,kz2
               jjz=jz+kz
               if(jjz>nz.or.jjz<=0) cycle
               do kx=-nkern,nkern
                  jjx=jx+kx
                  if(jjx>nx.or.jjx<=0) cycle
                  do ky=-nkern,nkern
                     jjy=jy+ky
                     if(jjy>ny.or.jjy<=0) cycle
                     s=s+img(jjx,jjy,jjz)*kern(kx,ky,kz)
                  end do
               end do
            end do
            tmpimg(jx,jy,jz)=real(s,4)
         end do
      end do
    end subroutine do_conv

  end subroutine conv_img
  
  !
  subroutine resort_sysmat
    integer ix,iy
    integer*2,allocatable,target::jxs(:),jys(:),jzs(:)
    integer*2,allocatable::jxse(:),jyse(:),jzse(:)
    integer*2,allocatable::jxso(:),jyso(:),jzso(:)
    real*4,allocatable::ws(:),wse(:),wso(:)
    integer*2,pointer::js(:)
    integer n,p,q,jo,je,k
    real s
    n=maxval(sysmat%n)
    allocate(jxs(n),jys(n),jzs(n),ws(n))
    allocate(jxse(n),jyse(n),jzse(n),wse(n))
    allocate(jxso(n),jyso(n),jzso(n),wso(n))
    do iy=1,my
       do ix=1,mx
          n=sysmat(ix,iy)%n
          if(n==0) cycle
          p=sysmat(ix,iy)%p
          if(p==0) cycle
          call mcps(loc(jxs),p,      n*2)
          call mcps(loc(jys),p+n*2,  n*2)
          call mcps(loc(jzs),p+n*2*2,n*2)
          call mcps(loc(ws), p+n*2*3,n*4)
          s=sum(ws(1:n))
          call free(p) ! to collect p and q
          !
          ! sorting along x
          !
          js=>jxs
          call separ
          sysmat(ix,iy)%nox=jo
          p=malloc(n*SIZE_OF_SYSMATEL)
          sysmat(ix,iy)%p=p
          call stob(p)
          !
          ! sorting along y
          !
          js=>jys
          call separ
          sysmat(ix,iy)%noy=jo
          q=malloc(n*SIZE_OF_SYSMATEL)
          sysmat(ix,iy)%q=q
          call stob(q)
       end do
    end do
    deallocate(jxso,jyso,jzso,wso)
    deallocate(jxse,jyse,jzse,wse)
    deallocate(jxs,jys,jzs,ws)
    !
  contains
    !
    subroutine separ
      jo=0
      je=0
      do k=1,n
         if(iand(js(k),1)==0) then
            je=je+1
            jxse(je)=jxs(k)
            jyse(je)=jys(k)
            jzse(je)=jzs(k)
            wse(je)=ws(k)
         else
            jo=jo+1
            jxso(jo)=jxs(k)
            jyso(jo)=jys(k)
            jzso(jo)=jzs(k)
            wso(jo)=ws(k)
         end if
      end do
    end subroutine separ
    !
    subroutine stob(dst)
      integer,intent(in)::dst
      if(jo>0) then
         jxs(1:jo)=jxso(1:jo)
         jys(1:jo)=jyso(1:jo)
         jzs(1:jo)=jzso(1:jo)
         ws (1:jo)= wso(1:jo)
      end if
      if(je>0) then
         jxs(jo+1:n)=jxse(1:je)
         jys(jo+1:n)=jyse(1:je)
         jzs(jo+1:n)=jzse(1:je)
         ws (jo+1:n)= wse(1:je)
      end if
      call mcps(dst,      loc(jxs),n*2)
      call mcps(dst+n*2,  loc(jys),n*2)
      call mcps(dst+n*2*2,loc(jzs),n*2)
      call mcps(dst+n*2*3,loc(ws), n*4)
    end subroutine stob
    !
  end subroutine resort_sysmat
  !
  subroutine set_idim(ifo,matrix_type) 
    use duty
    type(image_info),intent(inout)::ifo
    integer,intent(in)::matrix_type
    ifo%slice_width=0
    ifo%slice_separation=0
    ifo%number_of_slices=0
    ifo%width=0_2
    ifo%height=0_2
    ifo%dx=0
    ifo%dy=0
    ifo%bias=0
    ifo%intcpt=0  
    select case(matrix_type)
    case(MAT_IMG)
       ifo%slice_width=dz
       ifo%slice_separation=dz
       ifo%number_of_slices=nz
       ifo%width=int(nx,2)
       ifo%height=int(ny_true,2)
       ifo%dx=dv
       ifo%dy=dv
    case(MAT_PRJ)
       ifo%slice_width=cd      ! phony
       ifo%slice_separation=cd ! phony
       ifo%number_of_slices=mx*my
       ifo%width=int(mx,2)
       ifo%height=int(my,2)
       ifo%dx=cd
       ifo%dy=cd       
    end select
  end subroutine set_idim

  integer function mSiddong(dy2)
!USE OMP_LIB
    integer,intent(in)::dy2
    integer dx2
    real,pointer::X1(:)   ! start of LOR
    real,pointer::X2(:)   ! end of LOR
    real d12     ! length of LOR
    real Xp1(3)  ! coordinate of the first plane 
    real XpN(3)  ! coordinate of the last plane 
    ! a must be larger than 0 and less than 1.
    real wbuf(nbuf),wwbuf(nbuf) ! <<<
    integer jbuf(nbuf),per(nbuf)
    ! etc
    integer k,n
    integer,target::js(3)
    integer,pointer::jx,jy,jz
    integer ix2,iy2
    real norm
    real,target::X1s(3,mzs),X2s(3,mzs)
    real,target::w1s(mzs),w2s(mzs)
    integer izs1,izs2
    integer*2 ax
    real ds(3)
    pointer(pax,ax)
    real,parameter::wmin=1.e-9_dp
    real x0,y0,z0,a0
    !--------------------------------!
    real w,wtmp(128)
    !--------------------------------!
    integer jxmaxn,jymaxn,jxminn,jyminn
    real,parameter::neighbor_width=1.5
    real x,y
    logical fake
    integer in

    fake=(nbuf==0)
    mSiddong=0
    call init
    iy2=1+dy2
    main: do dx2=mx-1,dy2,-1
       if(max_mdxysqr>=0) then
          if(dx2**2+dy2**2>max_mdxysqr) cycle
       end if
       ix2=1+dx2
       call set_LOR_ends(ix2,iy2,X2s,w2s)
       k=0
       do izs2=1,mzs
          X2=>X2s(:,izs2)
          do izs1=1,mzs
             X1=>X1s(:,izs1)
             d12=sqrt(sum((X2-X1)**2))
             norm=w1s(izs1)*w2s(izs2)*(D/d12)**3
             do jz=1,nz
                z0=(jz-jz0)*dz
                a0=1.0_dp/2.0_dp-z0/D
                x0=X1(1)+a0*(X2(1)-X1(1))
                y0=X1(2)+a0*(X2(2)-X1(2))

                jxmaxn=int(nint(x0+neighbor_width)/dv+jx0)
                jxminn=int(nint(x0-neighbor_width)/dv+jx0)
                jymaxn=int(nint(y0+neighbor_width)/dv+jy0)
                jyminn=int(nint(y0-neighbor_width)/dv+jy0)
                
                in=0
                do jy=max(jyminn,1),min(jymaxn,ny)
                   y=(jy-jy0)*dv                   
                   do jx=max(jxminn,1),min(jxmaxn,nx)
                      x=(jx-jx0)*dv
                      w=exp(-((x-x0)**2+(y-y0)**2)/(2.0_dp))
                      in=in+1
                      wtmp(in)=w
                   end do
                end do
                wtmp(1:in)=wtmp(1:in)/sum(wtmp(1:in))*norm*cwz(jz)
                in=0
                do jy=max(jyminn,1),min(jymaxn,ny)
                   do jx=max(jxminn,1),min(jxmaxn,nx)
                      in=in+1
                      w=wtmp(in)
                      IF(STORE_JSW()==-1) EXIT MAIN
                   end do
                end do
             end do
          end do
       end do
       n=k
       if(fake) exit main
       if(n>0) then
          call add_a
          call alloc_sysmat_buf(n,sysmat(ix2,iy2))
          call unpack2(n,loc(jbuf),sysmat(ix2,iy2)%p,loc(wbuf))
       end if
    end do main
    
    if(fake) mSiddong=n*murphy_factor
    !
  contains 
    
    integer function store_jsw()
      store_jsw=0 ! OK
      if(jz<nzlo.or.jz>nzup &
           .or.jx<jxmin.or.jx>jxmax&
           .or.jy<jymin.or.jy>jymax) return
      if(w<=wmin) return
      k=k+1
      if(fake) return
      if(k>nbuf) then
         mSiddong=EID_BUFOVER_2
         store_jsw=-1
         return
         !     exit main
      end if
      pax=loc(jbuf(k))
      ax=transfer(jx,ax)
      pax=pax+2
      ax=transfer(jy,ax)
      pax=pax+2
      ax=transfer(jz,ax)
      wbuf(k)=w*sysmat_scale
    end function store_jsw
    !
    subroutine init
      Xp1(1)=-real(nx)*dv/2.0_dp
      Xp1(2)=-real(ny)*dv/2.0_dp
      Xp1(3)=-real(nz)*dz/2.0_dp
      XpN=-Xp1
      !
      jx=>js(1)
      jy=>js(2)
      jz=>js(3)
      !
      call set_LOR_ends(1,1,X1s,w1s)
      X1s(3,:)=-X1s(3,:)
      !
      ds(1:2)=dv
      ds(3)  =dz
    end subroutine init
    
    subroutine add_a
      integer src,trg,dst
      call persort_i(n,jbuf(1:n),per(1:n))
      src=1
      trg=1
      dst=1
      wwbuf(1:n)=wbuf(1:n)
      do while(src<n)
         src=src+1
         if(jbuf(src)/=jbuf(trg)) then
            dst=dst+1
            trg=src
            wwbuf(dst)=wbuf(src)
            jbuf(dst)=jbuf(src)
         else
            wwbuf(dst)=wwbuf(dst)+wbuf(per(src))
         end if
      end do
      n=dst
      wbuf(1:n)=wwbuf(1:n)
    end subroutine add_a
    !
  end function mSiddong

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine load_Ci(f,ci)
    character*(*),intent(in)::f
    real*4,intent(out)::ci(:,:,:,:)
    integer istat
    open(unit=1,file=f,iostat=istat,status="old",access="stream")
    if(istat/=0) call echo_red(mess_ioerr(EID_OPEN,f,istat))
    read(1,iostat=istat) ci
    if(istat/=0) call echo_red(mess_ioerr(EID_READ,f,istat))    
    close(1)
    call echo("loaded Ci: "//trim(f))
  end subroutine load_Ci

  subroutine save_Ci(f,v)
    character*(*),intent(in)::f
    real*4,intent(in)::v(:,:,:,:)
    character(MAXPATHLEN) fout
    integer istat
    fout=trim(f)//"_"//trim(i2a(size(v,1)))//"x"&
         //trim(i2a(size(v,2)))//".ci"
    call write_v(fout,-1,-1,loc(v),size(v),istat,OFMT_RAW)
    if(istat/=0) call echo_red(mess_ioerr(EID_WRITE,fout,istat))
  end subroutine save_Ci

  subroutine atcUProj(h)
    real*4,intent(inout)::h(:,:,:,:)
    real,parameter::mu=9.607e-3_dp ! mu of I=I_0*exp(-mu*x) in 1/mm
    integer ix1,iy1
    integer ddx,ddy,dx,dy
    real a,s
    s=sum(h)
    do ddy=1-my,my-1
       if(ddy>=0) then
          dy=ddy+1
       else
          dy=my+ddy+1
       end if
       do ddx=1-mx,mx-1
          if(ddx>=0) then
             dx=ddx+1
          else
             dx=mx+ddx+1
          end if
          a=exp(mu*sqrt(4*(ddx**2+ddy**2)+D**2)) ! <<<< `4' should be genlized 
          do iy1=1,my
             do ix1=1,mx
                h(ix1,iy1,dx,dy)=real(a*h(ix1,iy1,dx,dy),4)
             end do
          end do
       end do
    end do
    h=h*real(s/sum(h),4)
  end subroutine atcUProj
  
  !
  subroutine flipHistB1(h)
    real*4,intent(inout)::h(:,:,:,:)
    real*4 buf(mx,my,mx,my)
    integer ix1,iy1,ix2,iy2,im1,xoff,xoffr
    do iy2=1,my_true
       do ix2=1,mx
          do iy1=1,my
             do im1=1,md
                xoff=(mxb+mg)*(im1-1)
                xoffr=(mxb+mg)*((md-im1+1)-1)
                do ix1=xoff+1,xoff+mxb                            
                   buf(ix1,iy1,ix2,iy2)=h((ix1-xoff)+xoffr,iy1,ix2,iy2)
                end do
             end do
          end do
       end do
    end do
    h=buf
  end subroutine flipHistB1

  !
  subroutine flipHistX1(h)
    real*4,intent(inout)::h(:,:,:,:)
    real*4 buf(mx,my,mx,my)
    integer ix1,iy1,ix2,iy2,im1,xoff
    do iy2=1,my_true
       do ix2=1,mx
          do iy1=1,my
             do im1=1,md
                xoff=(mxb+mg)*(im1-1)
                do ix1=xoff+1,xoff+mxb                            
                   buf(ix1,iy1,ix2,iy2)=h(mxb-(ix1-1-xoff)+xoff,iy1,ix2,iy2)
                end do
             end do
          end do
       end do
    end do
    h=buf
  end subroutine flipHistX1
  !
  subroutine flipHistX2(h)
    real*4,intent(inout)::h(:,:,:,:)
    real*4 buf(mx,my,mx,my)
    integer ix1,iy1,ix2,iy2,im1,xoff
    do iy2=1,my_true
       do ix1=1,mx
          do iy1=1,my
             do im1=1,md
                xoff=(mxb+mg)*(im1-1)
                do ix2=xoff+1,xoff+mxb                            
                   buf(ix1,iy1,ix2,iy2)=h(ix1,iy1,mxb-(ix2-1-xoff)+xoff,iy2)
                end do
             end do
          end do
       end do
    end do
    h=buf
  end subroutine flipHistX2

  !
  subroutine flipHistY1(h)
    real*4,intent(inout)::h(:,:,:,:)
    real*4 buf(mx,my,mx,my)
    integer ix1,iy1,ix2,iy2    
    buf=0.0_4
    do iy2=1,my_true
       do ix2=1,mx             
          do iy1=1,my_true
             do ix1=1,mx
                buf(ix1,iy1,ix2,iy2)=h(ix1,my_true-iy1+1,ix2,iy2)
             end do
          end do
       end do
    end do
    h=buf
  end subroutine flipHistY1

end module PLANAR_PEM
